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Abstract 

We calculate the dynamic structure factor S(k,u) of liquid Ge (£-Ge) at 
temperature T = 1250 K, and of amorphous Ge (a-Ge) at T = 300 K, us- 
ing ab initio molecular dynamics. The electronic energy is computed using 
density-functional theory, primarily in the generalized gradient approxima- 
tion, together with a plane wave representation of the wave functions and 
ultra-soft pseudopotentials. We use a 64-atom cell with periodic boundary 
conditions, and calculate averages over runs of up to about 16 ps. The cal- 
culated liquid S(k,uj) agrees qualitatively with that obtained by Hosokawa 
et al, using inelastic X-ray scattering. In a-Ge, we find that the calculated 
S(k,u) is in qualitative agreement with that obtained experimentally by Ma- 

* Present address: Institute for Physical Science and Technology, University of Maryland, College 
Park, Maryland 20742; email jdchai@wam.umd.edu 
f Corresponding author: tel. (614)292-8140; fax (614)292-7557; email stroud® mps.ohio-state.edu 

1 



ley et al. Our results suggest that the ab initio approach is sufficient to allow 

approximate calculations of S(k, uS) in both liquid and amorphous materials. 
PACS numbers: 61.20.Gy, 61.20.Lc, 61.20.Ja, 61.25.Mv, 61.43.Dq 
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I. INTRODUCTION 



Ge is a well-known semiconductor in its solid phase, but becomes metallic in its liquid 
phase. Liquid Ge (£-Ge) has, near its melting point, an electrical conductivity characteristic 
of a reasonably good metal (~ 1.6 x lCT^^cm -1 [1]), but it retains some residual structural 
features of the solid semiconductor. For example, the static structure factor S(k), has a 
shoulder on the high-/c side of its first (principal) peak, which is believed to be due to a 
residual tetrahedral short-range order. This shoulder is absent in more conventional liquid 
metals such as Na or Al, which have a more close-packed structure in the liquid state and a 
shoulderless first peak in the structure factor. Similarly, the bond-angle distribution function 
just above melting is believed to have peaks at two angles, one near 60° and characteristic 
of close packing, and one near 108°, indicative of tetrahedral short range order. This latter 
peak rapidly disappears with increasing temperature in the liquid state. 

These striking properties of £-Ge have been studied theoretically by several groups. Their 
methods fall into two broad classes: empirical and first-principles. A typical empirical cal- 
culation is that of Yu et al [2], who calculate the structural properties of £-Ge assuming 
that the interatomic potentials in £-Ge are a sum of two-body and three-body potentials 
of the form proposed by Stillinger and Weber [3]. These authors find, in agreement with 
experiment, that there is a high-k shoulder on the first peak of S(k) just above melting, 
which fades away with increasing temperature. However, since in this model all the poten- 
tial energy is described by a sum of two-body and three-body interactions, the interatomic 
forces are probably stronger, and the ionic diffusion coefficient is correspondingly smaller, 
than their actual values. 

In the second approach, the electronic degrees of freedom are taken explicitly into ac- 
count. If the electron- ion interaction is sufficiently weak, it can be treated by linear response 
theory [4]. In linear response, the total energy in a given ionic configuration is a term which 
is independent of the ionic arrangement, plus a sum of two-body ion- ion effective interac- 
tions. These interactions typically do not give the bond-angle-dependent forces which are 
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present in the experiments, unless the calculations are carried to third order in the electron- 
ion pseudopotential [4]. Such interactions are, however, included in the so-called ab initio 
approach, in which the forces on the ions are calculated from first principles, using the 
Hellman-Feynman theorem together with density-functional theory [5] to treat the energy 
of the inhomogeneous electron gas. This approach not only correctly gives the bond-angle- 
dependent ion-ion interactions, but also, when combined with standard molecular dynamics 
techniques, provides a good account of the electronic properties and such dynamical ionic 
properties as the ionic self-diffusion coefficients. 

This combined approach, usually known as ab initio molecular dynamics, was pioneered 
by Car and Parrinello [6], and, in somewhat different form, has been applied to a wide range 
of liquid metals and alloys, including £-Ge [7-9], £-Ga x Gei_ x [10], stoichiometric III-V ma- 
terials such as £-GaAs, £-GaP, and £-lnP [11,12], nonstoichiometric ^-Ga^As]^ [13], £-CdTe 
[14], and £-ZnTe [15], among other materials which are semiconducting in their solid phases. 
It has been employed to calculate a wide range of properties of these materials, including 
the static structure factor, bond-angle distribution function, single-particle electronic den- 
sity of states, d. c. and a. c. electrical conductivity, and ionic self-diffusion coefficient. The 
calculations generally agree quite well with available experiment. 

A similar ab initio approach has also been applied extensively to a variety of amorphous 
semiconductors, usually obtained by quenching an equilibrated liquid state from the melt. 
For example, Car, Parrinello and their collaborators have used their own ab initio approach 
(based on treating the Fourier components of the electronic wave functions as fictitious 
classical variables) to obtain many structural and electronic properties of amorphous Si 
[16,17]. A similar approach has been used by Lee and Chang [18]. Kresse and Hafner [7] 
obtained both S(k) and g(r), as well as many electronic properties, of a-Ge, using an ab 
initio approach similar to that used here, in which the forces are obtained directly from the 
Hellmann-Feynman theorem and no use is made of fictitious dynamical variables for the 
electrons, as in the Car-Parrinello approach. A similar calculation for a-Si has been carried 
out by Cooper et al [19], also making use of a plane wave basis and treating the electron 
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density functional in the generalized gradient approximation (GGA) [20]. More recently, a 
number of calculations for a-Si and other amorphous semiconductors have been carried out 
by Sankey et al [21], and by Drabold and collaborators [22]. These calculations use ab initio 
molecular dynamics and electronic density functional theory, but in a localized basis. A 
recent study, in which S(k) and g(r) were computed for several ab initio structural models 
of a-Si, has been carried out by Alvarez et al [23]. 

Finally, we mention a third approach, intermediate between empirical and ab initio 
molecular dynamics, and generally known as tight-binding molecular dynamics. In this 
approach, the electronic part of the total energy is described using a general tight-binding 
Hamiltonian for the band electrons. The hopping matrix elements depend on separation 
between the ions, and additional terms are included to account for the various Coulomb 
energies in a consistent way. The parameters can be fitted to ab initio calculations, and 
forces on the ions can be derived from the separation-dependence of the hopping matrix 
elements. This approach has been used, e. g., to treat £-Si [24], a-Si [25], and liquid 
compound semiconductors such as £-GaAs and £-GaSb [26]. Results are in quite good 
agreement with experiment. 

In this paper, we extend the method of ab initio molecular dynamics to another dy- 
namical property of the ions: the dynamical structure factor, denoted S(k,u). While no 
fundamentally new theory is required to calculate S(k,u>), this quantity provides additional 
information about the time-dependent ionic response beyond what can be extracted from 
other quantities. The present work appears to be the first to calculate S(k, oo) using ab initio 
molecular dynamics. Here, we will calculate S(k, u) for £-Ge, where some recent experiments 
[27] provide data for comparison, and also for amorphous Ge (a-Ge). In the latter case, us- 
ing a series of approximations described below, we are able to infer the vibrational density 
of states of as-quenched a-Ge near temperature T = 300K in reasonable agreement with 
experiment. The calculated S(k,u) for the liquid also agrees quite well with experiment, 
especially considering the computational uncertainties inherent in an ab initio simulation 
with its necessarily small number of atoms and limited time intervals. 
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The remainder of this paper is organized as follows. A brief review of the calculational 
method is given in Section II. The results are presented in Section III, followed by a discussion 
and a summary of our conclusions in Section IV. 

II. METHOD 

Our method is similar to that described in several previous papers [9,10,13], but uses 
the Vienna Ab Initio Simulation Package (VASP), whose workings have been extensively 
described in the literature [28]. Briefly, the calculation involves two parts. First, for a given 
ionic configuration, the total electronic energy is calculated, using an approximate form of 
the Kohn-Sham free energy density functional, and the force on each ion is also calculated, 
using the Hellmann-Feynman theorem. Second, Newton's equations of motion are integrated 
numerically for the ions, using a suitable time step. The process is repeated for as many time 
steps as are needed to calculate the desired quantity. To hold the temperature constant, we 
use the canonical ensemble with the velocity rescaled at each time step. Further details of 
this approach are given in Ref. [9]. 

The VASP code uses ultrasoft Vanderbilt pseudopotentials [29], a plane wave basis for 
the wave functions, with the original Monkhorst-Pack (3x3x3) k-space meshes [30] and a 
total of 21,952 plane waves, corresponding to an energy cut-off of 104. 4eV. We use a finite- 
temperature version of the Kohn-Sham theory [31] in which the electron gas Helmholtz 
free energy is calculated on each time step. This version also broadens the one-electron 
energy levels to make the k-space sums converge more rapidly. Most of our calculations are 
done using the generalized gradient approximation (GGA) [20] for the exchange-correlation 
energy (we use the particular form of the GGA developed by Perdew and Wang [20]), but 
some are also carried out using the local-density approximation (LDA). 

In our iteration of Newton's Laws in liquid Ge (£-Ge), we typically start from the diamond 
structure (at the experimental liquid state density for the temperature of interest), then 
iterate for 901 time steps, each of 10 fs, using the LDA. To obtain S(k) within the GGA, 
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we start from the LDA configuration after 601 time steps, then iterate using the GGA for 
an additional 1641 10-fs time steps, or 16.41 ps. We calculate the GGA S(k) by averaging 
over an interval of 13.41 ps within this 16.41 time interval, starting a time t 2 after the start 
of the GGA simulation. We average over all t 2 , s from 1.0 to 3.0 ps. 

For comparison, we have also calculated S(k) within the LDA. This S(k) is obtained by 
averaging over 601 time steps of the 901 time-step LDA simulation. This 601-step interval 
is chosen to start a time t\ after the start of this simulation; the calculated LDA S(k) is also 
averaged over all t\S from lps to 3ps. 

To calculate quantities for amorphous Ge (a-Ge), we start with Ge in the diamond 
structure at T = 1600K and at the experimental density for that temperature, as quoted 
in Ref. [9]. Next, we quench this sample to 300 K, cooling at a uniform rate over, so as to 
reach 300 K in about 3.25 ps (in 10-fs time steps). Finally, starting from T = 300 K, we 
iterate for a further 897 time steps, each of 10 fs, or 8.97 ps, using the LDA. The LDA S(k) 
is then obtained by averaging over 5.97 of those 8.97 ps, starting a time t x after the system 
is reached 300K; we also average this S(k) over all t\% from 1.0 to 3.0 ps. To obtain S(k) 
within the GGA, we start the GGA after 5.7 ps of the LDA simulation, and then iterate 
using the GGA for an additional 18.11 ps in 10-fs time steps. The GGA S(k) is obtained 
by averaging over a 15.11 ps time interval of this 18.11 ps run, starting a time t 2 after the 
start of the GGA simulation; we also average over all values of t 2 from 1.0 to 3.0 ps. 

The reader may be concerned that the 3.25 ps quench time is very short, very unrepre- 
sentative of a realistic quench, and very likely to produce numerous defects in the quenched 
structure, which are not typical of the a-Ge studied in experiments. In defense, we note that 
some of these defects may be annealed out in the subsequent relaxation at low temperatures, 
which is carried out before the averages are taken. In addition, the static structure factor 
we obtain agrees well with experiment, and the dynamic structure factor is also consistent 
with experiment, as discussed below. Thus, the quench procedure appears to produce a 
material rather similar to some of those studied experimentally. Finally, we note that exper- 
iments on a-Ge themselves show some variation, depending on the exact method of sample 
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preparation. 

We have used the procedure outlined above to calculate various properties of £-Ge and a- 
Ge. Most of these calculated properties have been described in previous papers, using slightly 
different methods, and therefore will be discussed here only very briefly [32]. However, our 
results for the dynamic structure factor S(k, uj) as a function of wave vector k and frequency 
uj are new, and will be described in detail. We also present our calculated static structure 
factor S(k), which is needed in order to understand the dynamical results. 

S(k) is defined by the relation 

s ( k ) = ^<E ex P[* k • " ';(*))]>* - Wk,o, (!) 

where is the position of the i th ion at time t, N is the number of ions in the sample, and the 
triangular brackets denote an average over the sampling time. In all our calculations, we have 
used a cubic cell with N = 64 and periodic boundary conditions in all three directions. This 
choice of particle number and cell shape is compatible with any possible diamond-structure 
Ge within the computational cell. 

S(k, uj) is defined by the relation (for k ^ 0, w ^ 0) 

1 r°° 

S(k,u) = — y_ oo exp(^)(p(k,t)p(-k,0))rft, (2) 
where the Fourier component p(k, t) of the number density is defined by 

N 

p{k,t) = 5>xpHk-r,(t)). (3) 

i=i 

In our calculations, the average (...) is computed as 

(p(k, t)p(-k, 0)) = — J o p(k, h + t)p(-k, tjdh (4) 

over a suitable range of initial times t±. Because of the expected isotropy of the liquid or 
amorphous phase, which should hold in the limit of large N, S^k, a;) should be a function 
only of the magnitude k rather than the vector k, as should the structure factor S(k). 

Our calculations are carried out over relatively short times. To reduce statistical errors, 
we therefore first calculate 
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1 r'2 

S(k, u, ti, t 2 ) = — T7 / cftp(k, ti + t)p(-k, *i) exp(iw*). (5) 

7TiV Jti 

For large enough t 2 , S^k, a;, ti, t 2 ) should become independent of t 2 but will still retain some 
dependence on t±. Therefore, in the liquid, we obtain our calculated dynamic structure 
factor, S ca i c (k,tu), by averaging over a suitable range of ti from to Ati. 

1 /" A *i 

S calc (k,uj) = —— dt 1 S(k,u;,t 1 ,t2). (6) 

l\t\ Jo 

We choose our initial time in the t\ integral to be 1 ps after the start of the GGA calculation, 
and (in the liquid) Ati = 7 ps. We choose the final MD time t 2 = 16.41ps. For a-Ge, we 
use the same procedure but t 2 = 18. lips in our simulations. For our finite simulational 
sample, the calculated S(k) and S(k,u) will, in fact, depend on the direction as well as 
the magnitude of k. To suppress this finite-size effect, we average the calculated S(k) and 
S^k, a;) over all k values of the same length. This averaging considerably reduces statistical 
error in both S(k,u) and S(k). 

Finally, we have also incorporated the experimental resolution functions into our plotted 
values of S(k,cu). Specifically, we generally plot S av (k,u)/S(k), where S av (k,u) is obtained 
from the (already orient at ionally averaged) S(k, cu) using the formula 

/oo 
R(cu - cu')S(k,uj')duj', (7) 
-oo 

where the resolution function R(uj) (normalized so that Jf^, R(u)du = 1) is 

i?H = - ? i—exp(-^ 2 /cu 2 ). (8) 

In an isotropic liquid, we must have S(k,—cu) = S(k,cu), since our ions are assumed 
to move classically under the calculated first-principles force. Our orientational averaging 
procedure guarantees that this will be satisfied identically in our calculations, since for every 
k, we always include — k in the same average. We will nonetheless show results for negative 
ui for clarity, but they do not provide any additional information. 
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III. RESULTS 

A. S(k) for £-Ge and o-Ge 

In Fig. 1, we show the calculated S(k) for £-Ge at T = 1250 K, as obtained using 
the procedure described in Section II. The two calculated curves are obtained using the 
GGA and the LDA for the electronic energy-density functional; they lead to nearly identical 
results. The calculated S(k) shows the well-known characteristics already found in previous 
simulations [7,9]. Most notably, there is a shoulder on the high-/c side of the principal peak, 
which is believed to arise from residual short-range tetrahedral order persisting into the 
liquid phase just above melting. We also show the experimental results of Waseda et al 
[33]; agreement between simulation and experiment is good, and in particular the shoulder 
seen in experiment is also present in both calculated curves (as observed also in previous 
simulations). 

We have also calculated S(k) for a model of amorphous Ge (a-Ge) at 300 A". Our sample of 
a-Ge is prepared as follows. First, we create an equilibrated sample of £-Ge at a temperature 
of 1600A and the experimental density corresponding to that temperature (as quoted in 
Ref. [9]). Next, we quench this sample of £-Ge to 300A, by cooling at a uniform rate over 
a time interval of about 3.25ps. Finally, starting from T = 300A, we evaluate the structure 
factor by averaging over 800 time steps of 10 fs each, or 8 ps, after first discarding 100 time 
steps at T = 300K. We also average S(k) over different k vectors of the same length, as for 
£-Ge. We use the measured number density of 0.04370A~ 3 for a-Ge at T = 300A [34]. 

In Fig. 2, we show the calculated S(k) for a-Ge at T = 300, again using both the GGA 
and the LDA. The sample is prepared and the averages obtained as described in Section II. 
In contrast to £-Ge, but consistent with previous simulations [7,23], the principal peak in 
S(k) is strikingly split. The calculations are in excellent agreement with experiments carried 
out on as-quenched a-Ge at T = 300 K [34]; in particular, the split principal peak seen in 
experiment is accurately reproduced by the simulations. 
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We have also calculated a number of other quantities for both i-Ge and a-Ge, including 
pair distribution function g(r), and the electronic density of states n(E). 

For £-Ge, we calculated n(E) using the Monkhorst-Pack mesh with gamma point shifting 
(one of the meshes recommended in the VASP package). The resulting n(E) is generally 
similar to that found in previous calculations [7,9], provided that an average is taken over 
at least 5-10 liquid state configurations. Our n(E) for a-Ge [calculated using a shorter 
averaging time than that used below for S(k, u)] is also similar to that found previously [7]. 
Our calculated <?(r)'s for both i-Ge and a-Ge, as given by the VASP program, are similar 
to those found in Refs. [7] and [9]. The calculated number of nearest neighbors in the first 
shell is 4.18 for a-Ge measured to the first minimum after the principal peak in g(r). For 
i-Ge, if we count as "nearest neighbors" all those atoms within 3.4Aof the central atom (one 
of the cutoffs used in Ref. [7]) we find approximately 6.1 nearest neighbors, very close to 
the value of 5.9 obtained in Ref. [7] for that cutoff. Finally, we have recalculated the self- 
diffusion coefficient D(T) for i-Ge at T = 1250 K, from the time derivative of the calculated 
mean-square ionic displacement; we obtain a result very close to that of Ref. [9]. 

B. S(k,u) for i-Ge and a-Ge 

1. i-Ge 

Fig. 3 shows the calculated ratio S(k, uo)/ S(k) for i-Ge at T = 1250 K, as obtained using 
the averaging procedure described in Section II. We include a resolution function [eqs. (7) 
and (8)] of width hoo = 2.5 meV, the same as the quoted experimental width [27]. In Fig. 4, 
we show the same ratio, but without the resolution function (i. e., with ojq = 0). Obviously, 
there is much more statistical noise in this latter case, though the overall features can still 
be distinguished. 

To interpret these results, we first compare the calculated S(k, uo) in i-Ge with hydro- 
dynamic predictions, which should be appropriate at small k and cu. The This prediction 
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takes the form (see, for example, Ref. [35]): 

S(k,u) = 7-1 / 2D T k 2 \ 
71 S(k) 7 \uj 2 + (D T k 2 ) 2 J 

1 / Tk 2 Tk 2 

+ 7 {(uj + c s k) 2 + (Yk 2 ) 2 + (u - c s k) 2 + (rfc 2 ) 2 

Here 7 = cp/cy is the ratio of specific heats and constant pressure and constant volume, D T 

is the thermal diffusivity, c s is the adiabatic sound velocity, and Y is the sound attenuation 

constant. Dt and Y can in turn be expressed in terms of other quantities. For example, 

Dt = Kr/ipcp), where kt is the thermal conductivity and p is the atomic number density. 

Similarly, r = \ [0(7 — l)/7 + b], where a = nr/(pcv) and b is the kinematic longitudinal 

viscosity (see, for example, Ref. [35], pp. 264-66). 

Eq. (9) indicates that S(k, uj) in the hydrodynamic regime should have two propagating 
peaks centered at uj = ±c s k, and a diffusive peak centered at uj = and of width determined 
by D T . The calculated S(k, u)/S(k) for the three smallest values of k in Fig. 3, does show the 
propagating peaks. We estimate peak values of %oo ~ 10 meV for k = 5.60 nm" 1 , hu ~ 11 
meV for k = 7.92 nm -1 , and (somewhat less clearly) Twj ~ 13 meV for k = 9.70 nm -1 . The 
value of c s estimated from the lowest k value is c s ~ 2.7 x 10 5 cm/sec. (The largest of these 
three k values may already be outside the hydrodynamic, linear-dispersion regime.) 

These predictions agree reasonably well with the measured S(k, uj) obtained by Hosokawa 
et al [27], using inelastic X-ray scattering. For example, the measured sound-wave peaks 
for k = ±6 nm -1 occur near 10 meV, while those k = ±12 nm^ 1 occur at huo = 17.2 meV, 
Furthermore, the integrated relative strength of our calculated sound-wave peaks, compared 
to that of the central diffusion peak, decreases between k = 7.92 nm -1 and 12.5 nur 1 , 
consistent with both eq. (9) and the change in experimental behavior [27] between k — 6 
nm -1 and 12 nm -1 . 

Because S(k, uj) in Fig. 3 already includes a significant Gaussian smoothing function, a 
quantitatively accurate half-width for the central peak, and hence a reliable predicted value 
for Dt, cannot be extracted. A rough estimate can be made as follows. For the smallest k 
value of 5.6 nm -1 , the full width of the central peak at half-maximum is around 7.5 meV. If 
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the only broadening were due to this Gaussian smoothing, the full width would be around 
2hu = 5.5 meV. Thus, a rough estimate of the intrinsic full width is V7.5 2 — 5.5 2 = 5 
meV ~ 2hD T k 2 . This estimate seems reasonable from the raw data for S(k,cu) shown in 
Fig. 4. Using this estimate, one obtains D T ss 1.3 x 10~ 3 cm 2 /sec. 

The hydrodynamic expression for S(k,u)/S(k) was originally obtained without consid- 
eration of the electronic degrees of freedom. Since £-Ge is a reasonably good metal, one 
might ask if the various coefficients appearing eq. (9) should be the full coefficients, or just 
the ionic contribution to those coefficients. For example, should the value of D T which de- 
termines the central peak width be obtained from the full cp, cy, and k Ti or from only the 
ionic contributions to these quantities? For £-Ge, the question is most relevant for kt, since 
the dominant contribution to cp and cy should be the ionic parts, even in a liquid metal [4]. 
However, the principal contribution to k t is expected to be the electronic contribution. 

We have made an order-of-magnitude estimate of Dt using the experimental liquid num- 
ber density and the value Cp = (5/3) ks per ion, and obtaining the electronic contribution 
to Kp from the Wiedemann-Franz law [36] together with previously calculated estimates of 
the electronic contribution [9]. This procedure yields D T ~ 0.1 cm 2 /sec, about two orders 
of magnitude greater than that extracted from Fig. 3, and well outside the possible errors 
in that estimate. We conclude that the Dt which should be used in eq. (9) for £-Ge (and 
by inference other liquid metals) is the ionic contribution only. 

In support of this interpretation, we consider what one expects for S(k, u>) in a simple 
metal such as Na. In such a metal, ionic motions are quite accurately determined by effective 
pairwise screened ion-ion interactions [4]. Since the ionic motion is determined by such an 
interaction, the S(k, u) resulting from that motion should not involve the contribution of 
the electron gas to the thermal conductivity. Although £-Ge is not a simple metal, it seems 
plausible that its S(k, u) should be governed by similar effects, at least in the hydrodynamic 
regime. This plausibility argument is supported by our numerical results. 

For k beyond around 12 mn" 1 , the hydrodynamic model should start to break down, 
since the dimensionless parameter ojt (where r is the Maxwell viscoelastic relaxation time) 
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becomes comparable to unity. At these larger fc's, both our calculated and the measured [27] 
curves of S(k,u)/S(k) continue to exhibit similarities. Most notable is the existence of a 
single, rather narrow peak for k near the principal peak of S(k), followed by a reduction in 
height and broadening of this central peak as k is further increased. This narrowing was first 
predicted by de Gennes [37]. In our calculations, it shows up in the plot for k = 20.9 nm -1 , 
for which the half width of S(k,u)/S(k) is quite narrow, while at k — 28.5 and 30.7 nm -1 , 
the corresponding plots are somewhat broader and lower. By comparison, the measured 
central peak in S(k,u)/S(k) is narrow at k = 20 nm -1 and especially at k = 24 nm -1 , while 
it is broader and lower at k = 28 nm -1 [27]. 

The likely physics behind the de Gennes narrowing is straightforward. The half-width of 
S(k, uj) is inversely proportional to the lifetime of a density fluctuation of wave number k. 
If that k coincides with the principal peak in the structure factor, a density fluctuation will 
be in phase with the natural wavelength of the liquid structure, and should decay slowly, 
in comparison to density fluctuations at other wavelengths. This is indeed the behavior 
observed both in our simulations and in experiment. 

In further support of this picture, we may attempt to describe these fluctuations by a 
very oversimplified Langevin model. We suppose that the Fourier component p(k, t) [eq. 
(3)] is governed by a Langevin equation 

p(k,t) = -£p(k,t)+ V (t). (10) 

Here the dot is a time derivative, £ is a constant, and r](t) is a random time-dependent "force" 
which has ensemble average (rj) = and correlation function (rj(t)rj* (t')) — AS(t — t'). Eq. 
(10) can be solved by standard methods (see, e. g., Ref. [35] for a related example), with 
the result (for sufficiently large t) 

(p(k,t)p*(k,t + T)> = ^exp(-MO. (11) 

According to eq. (2), S(k, uj) is, to within a constant factor, the frequency Fourier transform 
of this expression, i. e. 
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/oo 
(irA/Z)exp(iuT)exp(-\T\£)dT, (12) 
-oo 

or, on carrying out the integral, 



7T A 

S(k,u)K^—. (13) 
oo z + 4 Z 



This is a Gaussian function centered at oo = and of half-width £. On the other hand, the 
static structure factor 

S(k) oc Lim r ^ (p(k, t)p*(k, t + r)) = —. (14) 

Thus, if the constant A is independent of k, the half-width £ of the function 5'(/c, u;) at wave 
number k is inversely proportional to the static structure S(k). This prediction is consistent 
with the "de Gennes narrowing" seen in our simulations and in experiment [27]. 

To summarize, there is overall a striking similarity in the shapes of the experimental and 
calculated curves for S(k,cu)/S(k) both in the hydrodynamic regime and at larger values of 
k. 



2. a-Ge 

We have also calculated the dynamic structure for our sample of a-Ge at T = 300K. 
The results for the ratio S(k,oo)/S(k) are shown in Figs. 5 and 6 for a range of k values, 
and, over a broader range of oo, in Fig. 7. Once again, both S(k, oo) and S(k) are averaged 
over different values of k of the same length, as described above. We have incorporated a 
resolution function of width Hojq = 2 meV into S(k,u). This width is a rough estimate for 
the experimental resolution function in the measurements of Maley et al [38]; we assume 
it to be smaller than the liquid case because the measured width of the central peak in 
S(k,u)/S(k) for a-Ge is quite small. 

Ideally, our calculated S(k, u) should be compared to the measured one. However, the 
published measured quantity is not S(k, oo) but is, instead, based on a modified dynamical 
structure factor, denoted G(k,oo), and related to S(k,oo) by [38] 
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c(^)=(p)( 



n(u,T) + 1 



) 



S(k, oj). 



(15) 



Here C is a k- and cj-independent constant, and n(oj,T) = \/[e* M l kBT — 1] is the phonon 
occupation number for phonons of energy E = hoj at temperature T. The quantity plotted 
by Maley et al [38] is an average of G(k,oj) over a range of k values from 40 to 70 nm _1 . 
These workers assume that this average is proportional to the vibrational density of states 
n v u,(u). The measured n V ib{oj) as obtained in this way [38] is shown in Fig. 8 for two 
different amorphous structures, corresponding to two different methods of preparation and 
having differing degrees of disorder. 

In order to compare our calculated S(k, oj) to experiment, we use eq. (15) to infer G(k, cu), 
then average over a suitable range of k. However, in using eq. (15), we use the classical form 
of the occupation factor, n(hoj) + 1 « ksT/hoj, This choice is justified because we have 
calculated S(k, oj) using classical equations of motion for the ions. We thus obtain for the 
calculated vibrational density of states 



In Fig. 8 we show two such calculated plots of n vi b(oj), as obtained by averaging eq. (16) 
over two separate groups of fc's, as indicated in the caption [39]. For comparison, we also 
show rimbioj) for a-Ge as calculated in Ref. [7] directly from the Fourier transform of the 
velocity- velocity autocorrelation function. 

The calculated plots for n V ib{oj) in Fig. 8 have some distinct structure, which arises from 
some corresponding high frequency structure in S(k,cu). The plot of n v n,{oj) for the group 
of smaller k's has two distinct peaks, near 8 meV and 29 meV, separated by a broad dip 
with a minimum near 18 meV. The plot corresponding to the group of larger k's has similar 
structure and width, but the dip is less pronounced. The two experimental plots also have 
two peaks separated by a clear dip. The two maxima are found around 10 and 35 meV, while 
the principal dip occurs near 16 meV. In addition, the overall width of the two densities of 
states is quite similar. 




(16) 
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The reasonable agreement between the calculated and measured n V ib{oj) suggests that our 
ab initio calculation of S(k, oj) for a-Ge is reasonably accurate. The noticeable differences 
probably arise from several factors. First, there are several approximations involved in going 
from the calculated and measured S(k,uys to the corresponding ra„ib(u;)'s, and these may 
be responsible for some of the discrepancies. Secondly, there may actually be differences 
between the particular amorphous structures studied in the experiments, and the quenched, 
then relaxed structure considered in the present calculations. (However, the similarities 
in the static structure factors suggest that these differences are not vast.) Finally, our 
calculations are carried out over relatively short times, using relatively few atoms; thus, 
finite-size and finite-time effects are likely to produce some additional errors. Considering 
all these factors, agreement between calculation and experiment is quite reasonable. 

Previous ab initio calculations for a-Ge [7] have also obtained a vibrational density of 
states, but this is computed directly from the ionic velocity-velocity autocorrelation function 
rather than from the procedure described here. The calculations in Ref. [7] do not require 
computing S(k,u). In the present work, by contrast, we start from S(k,u) (which is calcu- 
lated here for the first time in an ab initio calculation for a-Ge), and we work backwards to 
get n vi b(uj). In principle, our S(k, oj) includes all anharmonic effects on the vibrational spec- 
trum of a-Ge, though in extracting n v n,(uj) we assume that the lattice vibrates harmonically 
about the metastable atomic positions. In Fig. 8, we also show the results of Ref. [7] for 
n vi b(oj) as obtained from this correlation function. They are quite similar to those obtained 
in the present work, but have a somewhat deeper minimum between the two principal peaks. 

The quantity n vib {oj) could, of course, also be calculated directly from the force constant 
matrix, obtained by assuming that the quenched configuration is a local energy minimum 
and calculating the potential energy for small positional deviations from that minimum using 
ab initio molecular dynamics. This procedure has been followed for a-GeSe 2 , for example, by 
Cappelletti et al [40]. These workers have then obtained S(q,u>) versus q from their n v n,(uj) 
at selected values of oj, within a one-phonon approximation. However, as noted above, the 
present work produces the full S(k, oj) and thus has, in principle, more information than 
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n vib {to). 

IV. DISCUSSION AND CONCLUSIONS 

The present results show that ab initio molecular dynamics can be used to calculate the 
dynamic structure factor S(k, u>) for both liquid and amorphous semiconductors. Although 
the accuracy of the calculated S(k, u) is lower than that attained for static quantities, such 
as S(k), nonetheless it is sufficient for comparison to most experimental features. This is 
true even though our calculations are limited to 64-atom samples and fewer than 20 ps of 
elapsed real time. 

We have presented evidence that the calculated S(k, u)/S(k) in £-Ge agrees qualitatively 
with measured by inelastic X-ray scattering [27], and that the one calculated for a-Ge leads 
to a vibrational density of states qualitatively similar to the quoted experimental one [38]. 
Since such calculations are thus shown to be feasible, our work should spur further numer- 
ical studies, with longer runs on larger samples, to obtain even more detailed information. 
Furthermore, we can use these dynamical simulations to probe the underlying processes at 
the atomic scale which give rise to specific features in the measured and calculated S(k,u). 
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Figure Captions 

1. Static structure factor S(k) for l-Ge at T = 1250/T, just above the experimental 
melting temperature. Full curve: present work, as calculated using the generalized 
gradient approximation (GGA; see text). Dashed curve: present work, but using 
the local density approximation (LDA; see text). Open circles: measured S(k) near 
T = 1250 K, as given in Ref. [33]. 

2. Full curve: Calculated S(k) for a-Ge at T = 300K, as obtained using the GGA. 
Structure is prepared as described in the text. Open circles: measured S(k) for a-Ge 
at T = 3007C as given in Ref. [34]. 

3. Calculated ratio of dynamic structure factor S(k, uj) to static structure factor S{k) for 
l-Ge at T = 1250-K" for several values of k, plotted as a function of uj, calculated using 
ab initio molecular dynamics with a MD time step of 10 fs. For clarity, each curve 
has been vertically displaced by 0.05 units from the curve below. In each case, the 
plotted curve is obtained by averaging both the calculated S(k, u) and the calculated 
S(k) over all values of k of the same length. We also incorporate a Gaussian resolution 
function of half-width hu = 2.5 meV, as in eqs. (7) and (8). This value of u is chosen 
to equal the quoted experimental resolution [27]. 

4. Same as in Fig. 3, but without the resolution function. 

5. Calculated S(k,u)/S(k) for a-Ge at T = 300 K at k < 35 nn" 1 , plotted as a function 
of uj Again, each curve has been vertically displaced by appropriate amounts from 
the one below it, as evident from the Figure, and both S(k,uj) and S(k) have been 
plotted after an average over all k's of the same length. We also incorporate a Gaussian 
resolution function of half-width huj = 2 meV. This value is chosen to give the best 
results for n vi b(uj) as measured by Ref. [38]. The time step here is 10 fs. 

6. Same as Fig. 5, but at k > 35 nm _1 . 
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7. Calculated S(k,u)/S(k) as in Figs. 5 and 6 but including higher frequencies to. The 
vertical displacements are the same as in Figs. 5 and 6. At huo = 60 meV, the curves 
are arranged vertically in order of increasing frequency. 

8. Full curve: calculated vibrational density of states n v n,(ou), in units of 10 -3 states/meV. 
n vi b(u>) is obtained from the resolution-broadened S(k,cu) and S(k) of the previous 
Figure, using the formula n vi b(uj) = (G ca i c (k,ui)), where G ca i c (k,uj) is given by eq. 
(16), and the averaging is carried out over the three magnitudes of k near 40 nm" 1 for 
which we have computed S(k,cu). Dashed curve: same as full curve, but calculated 
by averaging over the six magnitudes of k near 90 nm -1 for which we have computed 
S(k,u). The open circles and open stars denote the measured n v n,(uj), as reported in 
Ref. [38] for two forms of a-Ge. Finally, the open diamonds denote n v n,(ui) as calculated 
in Ref. [7] from the ionic velocity- velocity autocorrelation function (dot-dashed curves). 
In all plots except that of Ref. [7], n V ib(ou) is normalized so that ^ max n(uj)d(huj) = 1. 
Umax is the frequency at which n vib {uj) — > 0, and is estimated from this Figure by 
extrapolating the right hand parts of the solid and dashed curves linearly to zero. 
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